Langevin theory of absorbing phase transitions with a conserved magnitude 
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The recently proposed Langevin equation, aimed to capture the relevant critical features of 
stochastic sandpiles, and other self-organizing systems is studied numerically. This equation is 
similar to the Reggeon field theory, describing generic systems with absorbing states, but it is cou- 
pled linearly to a second conserved and static (non-diffusive) field. It has been claimed to represent 
a new universality class, including different discrete models: the Manna as well as other sandpiles, 
reaction-diffusion systems, etc. In order to integrate the equation, and surpass the difficulties asso- 
ciated with its singular noise, we follow a numerical technique introduced by Dickman. Our results 
coincide remarkably well with those of discrete models claimed to belong to this universality class, in 
one, two, and three dimensions. This provides a strong backing for the Langevin theory of stochastic 
sandpiles, and to the very existence of this new, yet meagerly understood, universality class. 
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Aimed at shedding some light at the origin of "or- 
der'' in Nature, some different routes to organization 
have been proposed in the last fifteen years or so. In 
particular, the concept of self- organization, as exempli- 
fied by sandpiles 0, |2j 13 (for reviews see |3, IE Hi) 
(one of the canonical instances of self-organizing systems) 
has generated a rather remarkable outburst of interest. 
In order to rationalize sandpiles in particular, and selj- 
organized criticality (SOC) in general, and to understand 
their critical properties, it has been recently proposed 
to look at them as systems with many absorbing states 
11 0, II I1II3- The underlying idea is that in the ab- 
sence of external driving sandpile models get eventually 
trapped into stable configurations from which they can- 
not escape, i.e. absorbing states (AS) 0,0,0. In or- 
der to make the aforementioned connection more trans- 
parent, the notion of fixed energy sandpiles (FES) was 
introduced. These modified sandpiles are defined in such 
a way that they share the microscopic rules with their 
standard (slowly driven) counterparts, but with no driv- 
ing (no addition of sand-grains) nor dissipation; i. e. the 
total amount of sand (energy) becomes a conserved quan- 
tity that acts as a free tunable parameter. In this way, 
if a standard sandpile in its stationary, critical state has 
a density of grains (or energy) Cc, it can be verified that 
its fixed-energy counterpart exhibits a transition from an 
active to an absorbing phase at precisely Cc while it is 
in an absorbing (active) state below (above) that energy 
density. Slow driving and dissipation define a mechanism 
able to pin the system to its critical point 6, 7, 8, 9, l^]- 

Using this analogy with systems with AS, a field 



theoretical description of stochastic sandpiles has been 
proposedH 0,1113, which includes the two more rel- 
evant features of stochastic sandpiles: i) the presence 
of infinitely many AS and, ii) the global conservation 
of the total energy. The phcnomenological field theory 
(Langevin equation) aimed to capture the relevant criti- 
cal features of this type of systems is similar to the well 
known Reggeon field theory (RFT) ^12, 13, 14] (describ- 
ing generic systems with AS) but it is coupled linearly to 
a conserved energy field, namely 13 • 



dtp= Da^^p 



fip — Xp^ -\- ujpcj) -\- r]{x, t) 



(1) 
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where Da, Dc, fi, A and uj are constants, p{x,t) and 
4>{x, t) are the activity and the energy field respectively, 
7y is a zero-mean Gaussian noise with: 

{r^{x,t)r]{x' ,t')) = 2 p{x,t) S{x - x') S{t - t'). (2) 

As usual in systems with AS, the variance of 77 is propor- 
tional to p [ialT3|. 

Soon after the introduction of the previous Langevin 
equation its range of applicability was extended, as it was 
conjectured to describe all systems with many AS and 
an auxiliary conserved and non-diffusive (or static) field 
|15| . In particular for a reaction-diffusion model in this 
family an equation similar to Eq.(^ was derived rigor- 
ously by using standard Fock-space formalism techniques 
mm 113' To be more precise, we should mention that 
the derived set of equations, includes some higher order 
terms, as noise crossed-correlations, whose role in the 
asymptotic properties is unclear. 

A priori, it is not straightforward to decide from a 
field theoretical point of view whether the extra conser- 
vation law induces a critical behavior different from that 
of RFT or if, on the contrary, it is an irrelevant perturba- 
tion at the RFT renormalization group fixed point |l4j . 



2 



From the theory side, it has been recently argued by van 
Wijland that the Langevin equation is renormahzable in 
dc = 6 l2l] L while other authors have previously claimed 
dc = 4 (aJi, • Some mean field results and simulations 
in high dimensions of discrete models J^, and also a 
new method recently proposed by Lubeck to determine 
the upper critical dimension of systems with AS ^2d\ leads 
rather convincingly to dc = 4, but we are still far from a 
full clarification of these issues at a field theoretical level. 
In any case, it is commonly accepted, from numerical 
evidence, that this constitutes a new universality class, 
usually called Manna-class, or C-DP (in the spirit of Ho- 
henberg and Halperin |22)[i|i[ll|lllllllllll|23. 

In order to clarify the situation, it is the purpose of 
this paper to integrate numerically Eq.l^ in 1, 2, and 3 
dimensions. In this way we will verify whether this set of 
Langevin equations (without the higher order terms ob- 
tained in 'l6|) describes correctly the critical properties 
of the discrete models reported to belong to this class. 

With this purpose we employ an inte grat ion scheme in- 
troduced by Dickman some years back ^3l| which, to the 
best of our knowledge, is the only working method for 
Langevin equations including a RFT-like type of noise 
[i. e. with variance linear in the activity field). We 
will verify that indeed the Langevin equation as it is re- 
produces remarkably well the known exponents (as mea- 
sured in discrete models in this class), thus providing us 
with a sound base for further theoretical analyses of this 
universality class and of the role of conservation in self- 
organizing systems Q . 

The model- We integrate numerically Eq.(^ 
A technical problem appears when a standard (Euler 
|23 |) discretization scheme is used: due the symmetry 
of r] around zero, and the fact that the noise-term dom- 
inates the evolution whenever the density field is suffi- 
ciently small, negative (unphysical) local values of the 
density field can be generated. 

In order to overcome this difficulty a different, non- 
trivial, integration scheme was proposed by Dickman for 
the RET ^231 . It consists in discretizing the density field 
p as well as time and space. The quanta of density of 
activity can be taken proportional to the discrete time 
step, Ap = At, in such a way that the continuity of p is 
recovered in the thermodynamic limit {At, Ax —^ 0). The 
activity density at a given site i and time t is then given 
by p{i,t) = m{i,t)Ap, where m{i,t) takes integer values. 
Further details of the scheme can be found in 23] . It has 
been successfully applied to both the RET and to 
systems with many AS |2^, leading to good estimations 
of phase diagrams and critical properties. In order to 
extend the algorithm to our problem, the second equation 
of CQl is integrated using an usual Euler scheme with a 
continuous field (t){x,t). While for the equation of p, we 
follow a subtle different procedure based on Dickman's 
method. First, we calculate 

f{i, t + At)- f{i, t) = At[DaVlm{i, t) - pm{i, t) 
-XApm'^{i,t) +Lum{i,t),i){i,t)] + am^^^{i,t)'q' {i,t){3) 



where At = Ap, rj' is a zero-mean Gaussian white noise, 
is the discrete Laplacian operator and f(x,t) is an 
auxiliary continuous field. Then, after each integration 
step, the number of quanta of p, m{x,t), is updated ac- 
cording to: 

m{x, t + At) ^ m{x, t) 4- int[/(2;, t + At)] 

(4) 

f{x,t + At) f{x,t + At) - mt[ f{x,t + At)]. 

Initial conditions are taken as follows: i) (j){x, t — 0) = 
(/>o X (1 -f a\/'^e{x,t)), where e is a normahzed Gaus- 
sian noise with zero average, and a is a constant estab- 
lishing the range of relative variation allowed to (j) with 
respect to its mean value, (po- (l>o, is the control param- 
eter, and except for transient effects results should not 
depend on a. ii) The initial condition for p is chosen by 
randomly distributing active-field quanta, in such a way 
that p{x, t = 0) < 4>{x, t — 0) everywhere. 

We have carried out extensive simulations of the cou- 
pled Eqs.l^ in one, two and three dimensional lat- 
tices. In all the cases, the time-mesh has been fixed to 
At — 0.01, and Ax — 1, (values below which we have 
verified that results are not further affected). We also 
fix Da = Dc = 5 and p = X = uj = a = l. The noise 
amplitude a is taken different for the various dimensions 
in order to fix the transition in a reasonable (but arbi- 
trary) value of (f>Q: cr = 1 in ID, a = 0.5 in 2D, and 
a — 0.35 in iD. We have verified that the total energy 
is conserved within the considered precision, in all cases. 
The number of runs goes from 10^ up to 10^ depending 
on system size. 

Results. - As we vary 4>q, a continuous transition sepa- 
rating the absorbing from active phase, is observed at a 
critical threshold 0c. Below (above) this value the sys- 
tem is absorbing (active) states. The usual scaling laws: 

P (00 - ^c)", ? ~ (00 - (Pc)-"^, r ^ {cj^o- cj,c)-'«, 
where ^ (r) is the correlation length (time) are expected 
to hold 0, Ellil. This leads to the definition of the 
dynamic exponent as r ~ with z = v\\lv\^. It is also 
expected that at the critical point the density of activity 
presents a power law decay with time, p ^ t~^ . However, 
in some models, an anomalous critical time behavior of 
p has been reported [1 Elll E 113 ■ We shall later 
return to this issue. As usual the finite size of simulated 
systems induces the possibility of falling into the AS even 
for 00 > 0c. This fact has two consequences. The den- 
sity of activity p does not reach a stationary state close 
to the critical point. Hence, we are forced to consider the 
density of surviving trials, Psurv^ in order to realize the 
finite size analysis of p. On the other hand, this provides 
us with a method to measure the dynamic exponent z, by 
determining a characteristic decaying time as a function 
of system size. We have studied systems of linear size up 
to L = 4000 in ID, L = 400 in 2D and L = 80 in 31). The 
dependence of the stationary activity density on system 
size for several values of 0o in a one dimensional system 
is shown in Fig. 1. From this picture, we deduce the 
position of the critical point in ID, 0c(ll?) = 1.6371(1) 
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FIG. 1: Stationary value of psurv for different system sizes in 
ID. Squares correspond to (^0 = 1.6369, circles to 00 = 1.6371 
and diamonds to (j)o = 1.6373. In the inset, the same graph is 
displayed but for 2D data; squares are for <^o = 0.631, circles 
for — 0.6325 and diamonds for (j>o = 0.635. 
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FIG. 2; Evolution of p (continuous line) and paurv (dashed 
lines) for several system sizes in ID. The curve of p is for 
L — 4000 and those of psurv for from top to bottom L = 
20, 100 and 500. The slope of the straight line is d = 0.14. In 
the inset, the same data in 2D are represented; psurv curves 
correspond from top to bottom to L = 10, 25, 70 and L — 280 
and p to L = 280. The slope of the line is 6^ = -0.65 l27l| . 



(numbers in parentheses correspond to the statistic un- 
certainty in the last digit). Also, from the slope of the 
log-log plot, we obtain P/v±{lD) = 0.213(6). The ex- 
ponent /3 may be estimated in an independent way from 
the scaling of the stationary value of psurv for large sys- 
tem sizes, as a function of (^o — 0c) above the critical 
point. This gives P{1D) = 0.28(1). By studying the 
time evolution of the characteristic time of the surviving 
probability P{t) at criticality, we obtain z{lD) — 1.47(2). 
Finally, the exponent 0{lD) = 0.14(1) may be measured 
from the critical power law decay of p in time, as may be 
seen in Fig. 2. Errors in these exponents mostly come 
from the uncertainty in the determination of the criti- 
cal point. Repeating this process in higher dimensions. 



FIG. 3: Anomalous time decay of the activity density in 2D, 
for L = 280 and 0o = 0.711. 



we find M'^D) = 0.6325(5) and M^D) = 0.456(1), to- 
gether with the critical exponents listed in Table I. In 
the table, we have also included the critical exponents 
of discrete models claimed to belong to this same uni- 
versality, and also (for comparison) those of the Directed 
Percolation (DP) class. Observe the rather remarkable 
agreement (within errorbars) between all the measured 
exponents are their counterparts in discrete models. Let 
us remark that, for those exponents for which the differ- 
ences with DP- values are larger, our values also deviate 
from DP. As we have already mentioned, some models in 
the Manna universality class may present an anomalous 
behavior in the time decay of the activity density at the 
critical point H 113, El 111 II3| p{t, (t>c)- This anomaly 
implies th at, app arently, the scaling relation /? — Ov^^ 
fails Jj, 0,0, El 01 J and that p{t) may decay in a non- 
monotonous way at criticality. In our case, there is no 
anomalous decay in d = 1 (Fig. 2). However, the anomaly 
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present both in d 


= 2 and 
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As can 


be seen 
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Model 


/3 




2 




e 




Eq.CJ 


0.28(1) 


0.213(6) 


1.47(2) 


1.63(1) 


0.14(1) 


1 


C-DP 


0.29(2) 


0.217(9) 


1.55(3) 




0.140(5) 




DP 


0.276... 


0.252... 


1.580... 


1.733... 


0.159... 




Eq.lO 


0.66(1) 


0.85(8) 


1.51(3) 


1.27(7) 


0.50(5) 


2 


C-DP 


0.64(2) 


0.78(2) 


1.55(3) 


1.29(8) 


0.51(1) 




DP 


0.583(4) 


0.795(6) 


1.766(2) 


1.295(6) 


0.450(2) 




Eq.im 


0.84(5) 


1.44(5) 


1.69(3) 


1.07(8) 


0.93(3) 


3 


C-DP 


0.88(2) 


1.39(4) 


1.73(5) 


1.12(8) 


0.88(2) 




DP 


0.80(2) 


1.39(1) 


1.901(5) 


1.105(5) 


0.730(4) 



TABLE I: Critical exponents for steady state experiments in 
d = 1,2 and 3. Figures in parenthesis indicate the statisti- 
cal uncertainty in the last digit. C-DP exponents are from 
Refs. m [H] and DP exponents from Refs. In ID, 

13/ u± for C-DP has been calculated using the scaling relation 
f3/u± ^ zd. 
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in Fig. 3, the activity density decays initially faster than 
a power law, showing afterwards a non-monotonous be- 
havior. The fast decay explains the anomalous values of 
0, not satisfying scaling relations, usually reported in the 
literature j^]. Later on, p increases before reaching the 
steady state value, pstat{L) after a certain time tx{L) 
(calculated by fixing some arbitrary criterion). If the 
points (tx{L), pstatiL)) are represented in a log-log plot 
at the critical point, an alternative value for the expo- 
nent 9 is found, which is related to the saturation time 
scale. This value in two dimensional space, « 0.50 
(with a large statistic uncertainty), |23| is much closer to 
the more accurate measurements reported in the litera- 
ture for models in this class (6* = 0.51 for C-DP pl|V 
The common presence of anomalous behavior in discrete 
systems ||2i,iil^il^iiSa.iZ] and in the continuous theory re- 
inforces the claim that both belong to the same universal- 
ity class: they share not only the critical behavior but also 
the dynamical anomalies. A deeper study of the physical 
origin of this anomaly is still missing but, essentially, it is 
related to the existence of different time scales. 

Finally, let us mention that spreading experiments 
(performed using localized seeds of activity [ij, [13 ) ex- 
hibit some numerical instabilities, due to the presence of 
large gradients, that merit a separate analysis and are, 
therefore, left aside for a future work. 

Conclusions.- Strong numerical evidence confirms that 
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